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Abstract We consider a bilayer setup with two parallel planes of cold fermionic polar molecules when 
the dipole moments are oriented perpendicular to the planes. The binding energy of two-body states 
with one polar molecule in each layer is determined and compared to various analytic approximation 
schemes in both coordinate- and momentum-space. The effective interaction of two bound dimers is 
obtained by integrating out the internal dimer bound state wave function and its robustness under 
analytical approximations is studied. Furthermore, we consider the effect of the background of other 
fermions on the dimer state through Pauli blocking, and discuss implications for the zero-temperature 
many-body phase diagram of this experimentally realizable system. 



1 Introduction 

The cooling and trapping of ultracold polar molecules is a rapidly growing field at the moment [TJIHISl 
miSlinilZ] ) which generates huge theoretical interest [HUH] . To avoid heavy losses from chemical reactions 
[7] or many-body collapse 10 , trapping in low-dimensional geometries has been explored in multiple 
works [TT1[T^[T51IT3 . 15, 16 , 17, 18,19 , 20 . Very recently, a multilayer system of fermionic polar molecules 
was experimentally realized and chemical reactions rates that depend on the geometry and population 
of optical lattice states have been measured [5T]. In addition to the potential for ultracold controlled 
chemistry of such a system, the long-range character of the dipole-dipole interaction means that the 
interactions of the molecules can potentially support a number of exotic few- and many-body phases 

In the case of dipoles aligned perpendicular to the layers by an externally applied field the inter- 
action of dipoles in different layers has a very interesting structure with short-range attraction and 
long-range repulsion. In fact, in momentum space the interaction becomes negative definite in the two- 
dimensional (2D) limit where the layers are assumed infinitely thin [23] • This suggests that the system 
could exhibit pairing. However, the intralayer interaction for the same alignment is purely repulsive 
[42j . and this counteracts the pairing mechanism. This interplay of repulsive and attractive interac- 
tions in parallel planes is reminiscent of features in high-temperature superconductors and layers of 
cold polar molecules are therefore a very interesting model system. 

Here we focus on the case of a bilayer which already contains many of the features expected of the 
general many-layer system. The two-dimensional two-component Fermi gas with attractive short-range 
interactions is a well-studied problem. In particular, the interesting physics of the BCS-BEC crossover 
from weakly paired atoms to strongly bound dimers was obtained some years ago [i51BilH5] . It was 
shown that the crossover takes place when the Fermi energy of the system becomes comparable to the 
binding energy of the two-body bound state between the two components (typically electron spins) of 
the generic Fermi gas. The existence of the latter state is guaranteed by the famous Landau criterion 
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which states that a potential with a negative definite volume integral always supports a bound state 
[46] . A bilayer with fermionic polar molecules is similar to the 2D Fermi gas when mapping the spin to 
a layer index. The important role of the two-body bound state in the many-body physics of the bilayer 
is now clear. 

The bilayer setup with perpendicular dipole moments has a two-body bound state which is sup- 
ported by the interlayer potential whose inner part is the only source of attraction. However, the spatial 
integral of the potential turn out to be zero and Landau's criterion cannot be applied to secure a bound 
state at any coupling strength of the dipolar molecules. An early proof of the existence of a bound state 
even for vanishing spatial integral was given by Simon ^47) . In the initial stages of the theoretical work 
on dipolar systems in 2D, this result appears to have been forgotten and Gaussian approximations that 
predict a critical coupling strength for a bound two-body state to arise were used [iSlfT^I^I^ . The 
energy of the bound state for the interlayer potential of a bilayer with perpendicular dipoles was first 
calculated within a scattering approach in and by solving the Schrodinger equation in [5U1IM] . The 
highly non-trivial dependence of the bound state energy on the coupling strength was subsequently 
calculated for perpendicular [551f5Hl[51] , and for arbitrary orientation of the dipoles with respect to the 
planes j57] . 

In this paper we consider the properties of the two-body bound dimer for dipolar fermions in a 
bilayer with dipoles oriented perpendicular to the layers. The exact binding is obtained by numerical 
solution of the Schrodinger equation and different approximation schemes based on Gaussian and 
exponential wave functions are tested against the exact result. The momentum-space wave function is 
then calculated and used to discuss an effective dimer-dimer interaction for use in the corresponding 
many-body problem. Here we find that appropriate Gaussian wave functions provide good analytic 
approximations to the exact result everywhere except the weak-coupling limit where the wave function 
becomes very delocalized in space. Lastly, we consider the bound state in the context of the BCS-BEC 
crossover. We include the Pauli blocking effect of the fermionic background by solving the momentum- 
space Schrodinger equation and compare to a Gaussian approximation scheme. The binding energy of 
the dimer is found to decrease rapidly with increasing size of the background Fermi sea. The latter 
means that the (quasi)-BEC regime where the system effectively consists of bound dimers occupies 
only a small region of the phase diagram at very low density, making it hard to reach experimentally. 



2 Model 

We consider the case of a bilayer geometry with a single polar molecule of mass M and dipole moment 
d in each layer. We assume the molecules are fermions, but this will only play a role later when 
considering the effect of Pauli blocking on the binding energy. We assume that the dipole moments 
have been aligned by an external field so that they are oriented perpendicular to the layers in which 
the molecules move. In this setup the dipolar potential is 



V{r) = 



(r2 + d2)5/2 



(1) 



where d is the interlayer distance and = d? /Atteq (assuming electric dipole moments, d). The three- 



dimensional distance between the molecules is 



d , i.e. r is the relative distance in the plane. For 



r/d <C 1 the potential behaves as an oscillator, whereas for r/d 3> 1 it has a repulsive tail. Note 
the interesting feature that j d^rV{r) — 0, which implies a highly non-trivial behavior of the binding 
energy for small dipole moments 47. 50 , 55 . 541 (M1[571I5^ . 

Since the potential has cylindrical symmetry, the two-body Schrodinger equation can be decom- 
posed into partial waves !Z'(r) = ^m{T)'4'm{4') with radial equation 



2/i 



1 - 4m2 



Umir) + V{r)um{r) = Eumir), 



(2) 



where fi = M/2 is the reduced mass and Um{f) — y^r)Rjn{r) is the reduced radial wave function. In 
this paper we will only be concerned with the lowest m = state which yields the radial equation 



d^ 
dx^ 



1 jj ^ -'^ 
4a;2 (a;2 



1)5/2 



E 



uo{x) = 0, 



(3) 



Fig. 1 (color online) Schematic picture of the bilayer setup with one dipolar molecule in each layer at distance 
r. The dipole moment is polarized perpendicular to the bilayer planes by an externally applied field. 



where we have rescaled the equation to dimensionless units, i.e. x = r/d and E = Md?E/h^. The 
dimensionless strength of the dipolar interaction in the bilayer setup \sU = MD^/fi^d. Unless otherwise 
stated, we will use It? /Md^ as the unit of energy and d for lengths. 

The two-body Schrodinger equation in Eq. ([3| can be solved numerically by standard methods. 
We denote the exact energy obtained from solving this equation by Eq and the wave function by !Po- 
However, when the dimer appears as input for more involved calculations on the many-body physics of 
layered systems, it is extremely useful to have accurate analytic or semi-analytic approximations to the 
full wave function. This will be our main concern in this paper, and we now outline four approximation 
schemes that we want to compare. The first three are based on an increasingly more sophisticated 
oscillator approximation, whereas the last one is an exponential ansatz for the wave function. What 
will be evident is that it is very important for any approximation to properly reproduce the spatial 
profile of the wave function. 



2.1 Scheme 1 

The simplest and also crudest approximation is to expand the potential to quadratic order for r d, 

nr)^-^ + -^^r\ (4) 

which yields the relation iiu? = for the angular frequency of the oscillator. The corresponding 

wave function is simply 

1 jf_ 

where 1? — h/ jiuj. As we will see below, this gives a bad estimate for the binding energy, but does exhibit 
some of the right behavior for large U where the wave function localizes and becomes increasingly 
Gaussian. We denote the energy and wave function obtained from this first approximation scheme 
by El and tf^i. Recall that for Gaussian wave functions with parameter b in two dimensions we have 



2.2 Scheme 2 

Another strategy is to use the nice Gaussian wave function but employ the oscillator length corre- 
sponding to the correct binding energy. The relation between binding energy, Eb , and oscillator length 
parameter, I, is P = 2h^ /MEb- We now take Eb = \Eq\, i.e. we use the exact result for the energy 
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to determine the length scale of the Gaussian ground-state wave function '^2- This has the advantage 
that the wave function will become increasingly extended for small U and, in turn, small |-Eo|- 



2.3 Scheme 3 

The third use of Gaussian approximation schemes is somewhat more involved and builds on ideas 
discussed earlier in Ref. [SD] . We approximate the true potential by a model potential which is quadratic 
but with a constant (negative) shift. We then fit the oscillator frequency and the shift to ensure that 
(i) the binding energy is equal to the exact result and (ii) the node of the potential occurs at r/rf = 
which is the same position as for the real potential in Eq. ([T]). The latter condition reflects the idea that 
the inner attractive part is the important one that determines the properties of the system, whereas 
the outer repulsive tail of the true dipolar interaction has only limited influence. Assuming a potential 
of the form V{r) = Muj'^r^ /A — Eg, the two conditions yield the following matching relations between 
angular frequency, uj, shift. Eg > 0, and Eg 

]^M(fCj'^ = Es and Eq = HiL - Es- (6) 
These second order equations have two solutions 



x± = l±^l-2Eo, (7) 

where x± = M(Puj±/h and Eq — M(PEo/h^. The solutions have the properties that ^ 2 — Eq and 
X- Eq for \Eq\ — t' 0. We see that the x+ solution goes to a constant oscillator frequency, whereas 
the a;_ goes to a vanishing frequency. The approximating potential becomes 



V{r) ^ (x± jr^- 



1 

2 \d 



(8) 



where since Eq < the first factor is always positive and tends to zero for ji^ol — > in the weak- 
coupling limit. The solution a;„ is negative, and therefore = hx-/Md'^ < 0. Therefore this quantity 
cannot be interpreted as an oscillator frequency as it stands. However, if we take w_ = h\x-\/Md'^ , we 
get a positive frequency which has the desirable feature that the wave function is allowed to extend 
in space when \Eq\ is small. The price we pay is that now the binding energy in the oscillator is not 
equal to the exact result, Eq, but becomes 



Md^Q^/h - Md^Es/h^ ^2yi~~2Eo + Eo-2, (9) 

in dimensionless form. The energy is the uj- oscillator is then positive for —Eq < 4, but still goes 
to zero as \Eq\ — >■ 0. Since it is highly desirable to have an extended potential (i.e. a)_ — > 0) as 
|i?o| — we retain this solution despite the difference in ground state energy. The length scale of the 
Gaussian ground-state of the ui± oscillators has the simple form P/d'^ = 2/\x±\, which clearly diverges 
as |£'o| — ?> for the Cj- solution. We denote the corresponding wave function !p3. 

Below we will see that whereas the approximation using w_ compares poorly with exact results, 
using w+ gives a much better approximation in spite of its localized behaviour when approaching zero 
energy 2). The approximation with w+ has been used to calculate the bound state energy of 

chains with one molecule in each layer in the three- and four-layer cases, and yields results that are in 
excellent agreement with exact numerical methods [CTf5^[55] . 



2.4 Scheme 4 



The final approximation scheme uses an exponential wave function instead of a Gaussian but still uses 
the exact binding energy. This means that we reproduce the correct large-distance behavior of the 
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exact wave function. The wavenumber k, of the exponential is fixed similar to approximation scheme 
2, i.e. \Eq\ = Yi^K^ or k — y^M\Eo\/h?. The properly normalized wave function for this scheme is 

Mr) = (10) 

For this wave function we have (r^) — 3/2k^. Scheme 4 is special in the sense that it is designed to 
approximate the tail behavior of the full wave function. Since the exponential form is not correct at 
small distances it does not make sense to calculate the energy for this wave function. The correct wave 
function for bound states outside the range of the potential in two dimensions is the modified Bessel 
function of the second kind Kq, which does have the exponential tail but is logarithmically divergent 
at the origin. 



3 Dimer Properties 

We now compare the properties of the exact dimer solution from the Schrodinger equation to the 
various approximation schemes introduced above. We do this for the energy obtained in scheme 1, 
and for the wave function and radius squared for all schemes. Furthermore, we consider the Fourier 
transform of the wave functions and of their square which will be of use when discussing the effective 
interaction of dimers below. 



3.1 Binding Energy 

The first property of the dimer system to study is the binding energy. Since the ground-state has 
cylindrical symmetry, the exact solution can be obtained by integrating the Schrodinger equation and 
matching to a bound state wave function at large distances. This energy has been calculated by several 
groups in recent years [49 l [50 l [55 y 34| , [57] . In Fig. [2] the energy, Eq, is plotted as a function of the dipolar 
strength U. At very small t/ ~ 1 one sees the rapid decrease of the binding energy which to leading 
order goes as e~^^'^ [37[f551l54ll57lf5S] . The calculation was done for U > 0.9 as smaller U values have 
binding energies that were below the numerical precision. Note that the exact energy is recovered by 
definition in approximation schemes 2 and 3 above (i?2 = = Eq). 

Also shown in Fig. [2] is the result of calculating the energy by a crude harmonic approximation to 
the potential around r = 0. The energy has the following analytical expression 

^^Ei^V24U-2U. (11) 

For [/ < 6 we find Ei > 0. This is of course due to the zero-point energy which is very large in 
this particular approximation. However, at larger U we do see the correct slope of Ei as compared 
to the exact result, and the difference seems to be only an overall offset. This clearly indicates that 
approximating the wave function by a Gaussian should be an accurate description for larger U as noted 
in [50]. It is of course clear that this tight harmonic approximation will fail at small U. To remedy 
this situation, we now discuss a number of alternative approximations to the wave function which use 
knowledge of the exact result, but where the scale of the wave function varies with [/ in a different 
way as compared to the crude harmonic approximation of Wi . 

First, we check the energies calculated within the various approximation that are based on Gaussians 
where a nice analytical expression can be given for the total energy. In fact, we have 

^{H)=x'[l-Ufix)], (12) 

f{x) = -4a;2 + 2V^Erfc(a:)e"^' {2x^ + x), (13) 

where Erfc(a;) is the complementary error function. We used the quantity x — d/b, where b is the length 
parameter in the Gaussian wave function. The optimal variational energy is now found by finding the 
minimum of Eq. [12] with respect to x. This value can be compared to the choices of length parameter 
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Fig. 2 (color online) Dimer energy as a function of U in units of Mcf /fP for U > 0.9. The solid (black) line 
is the exact solution (-Bo), whereas the dashed (blue) line is the energy obtained by expanding the dipolar 
potential around the origin to second order in approximation scheme 1 (-Ei). The horizontal dotted (black) line 
indicates zero energy. 

based on the exact energy and the shape of the potential as outhned in approximation schemes 2 and 
3. We have chosen not to hst this approximation among the schemes above since it is not based on the 
knowledge of the binding energy for construction of the wave function which is what the comparison 
done here is about. The variational approach provides both a binding energy and a corresponding 
(Gaussian) wave function, which is as we will now see close to some of the other schemes. 

In Fig.|3]we plot the energies of the different Gaussian approaches along with the exact energy and 
also the energy result of the crude approximation in scheme 1. Note that all approximate solutions have 
positive energy below U 1.5. Comparing the solid (red) and dashed (black) curves we see that the 
variationally optimized solution is better than scheme 2 but not by a large margin, and the difference 
becomes small for large U where they both approach the exact result. The two dotted (green) curves 
from scheme 3 are interesting as w_ seems to do consistently worse than with the exception of 
small U where the a)+ solution goes positive. Note how scheme 2 and scheme 3 with ui+ does well in 
the large U and small U regimes respectively, where they agree also with the variational calculation. 



3.2 Wave function and Size 

Next we look at the wave function of the exact solution in comparison to the different approximation 
schemes. In Fig. |4]we plot all the different choices for strengths U = 1, 2, 5, and 10. As expected, 
scheme 1 is poor for small U and only becomes acceptable around U = 10, while the exponential 
approximation to the wave function is quite good in all cases. The former scheme simply can not 
deal with a rapid spreading of the wave function at small U. For schemes 2 and 3 we see quite good 
agreement with scheme 2 always giving a slightly better approximation. 

From the wave function we can determine the extension of the dimers in the various approximation 
by computing (r^). The results are shown in Fig. p^on a logarithmic scale which is convenient since 
at [/ ^ 1 the dimer becomes extremely extended. Again we see that the exponential form in scheme 4 
does remarkably well, whereas scheme 2 and 3 overestimate the size for U > 2.5. Scheme 1 has a very 
small radius for all U and is quite far off the exact result as was clear from the wave functions in Fig. [4j 
It is worth noting that all curves in Fig. [5] have similar decreasing behavior for large U although the 
asymptotic values are very different. It is also interesting that even though the exact result predicts a 
dimer of size less than the attractive part of the potential for C/ > 4, scheme 1 does not provide a good 
size estimate until much larger U even though it was built on only this inner attractive harmonic part 
of the full dipole-dipole potential. The lowest order expansion is therefore clearly insufficient. 
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Fig. 3 (color online) Dimer energy as a function of U in units of MS jf? for U > 0.9. The thick solid (blue) 
line is the exact solution (Eo) and the upper thin solid (black) line is the energy obtained by expanding the 
dipolar potential around the origin to second order in approximation scheme 1 (Ei). The solid (red) line is 
the energy based on the choice of length parameter in scheme 2, while the dashed (black) line is the optimum 
variational choice when using a Gaussian wave function. The two dotted (green) lines are for scheme 3 with the 
lower curve at large U corresponding to aj+ and the upper one to £D_ . The vertical dotted (black) line indicates 
zero energy. 

From the results above, we can conclude that the key issue to obtain a good approximation seems 
to be the choice of bound-state size. Schemes 2 and 4 which includes only the size estimated through 
the binding energy are always the more accurate choices, with scheme 3 which includes the zero of the 
original potential also doing well. The naive harmonic approximation to the potential seems to fail in 
all aspects as seen in Figures [3j |4] and [5j 



3.3 Fourier Transforms 

Below we are concerned with the effective interaction between dimers which can be calculated from the 
Fourier transform of the dimer wave function. We define the Fourier transform to momentum space as 



Since the wave function has cylindrical symmetry we can perform the angular integration and obtain 




(14) 




(15) 



where Jo{z) is the Bessel function of order zero. For the Gaussian wave function, 'P{x) = e 
the result is 



.2 



(16) 



while the exponential, !f'(x) 




yields 




(17) 



The remarks made about the various approximation schemes in coordinate space in relation to Fig. [4] 
hold for the Fourier transforms as well with scheme 1 doing quite bad while scheme 2,3, and 4 comparing 
much better with the exact result. 
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Fig. 4 (color online) Dimer radial wave function for (7 = 1, 2, 5 and 10. in units of MS" jTi" . The solid thick 
(blue) line is the exact solution. The other lines are approximation scheme 1 (solid red), 2 (dashed black), 3 
(dotted green), and 4 (dash-dotted brown). Note the difference in scales on the vertical axis. The upper (at 
r = 0) dotted green corresponds to 11)+ in scheme 3 and the lower one to £)_ . 
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Fig. 5 (color online) Mean radius squared, (r^), as function oil] on a logarithmic scale. The lines are as in 
Fig. li] The upper (at (7 = 1) dotted (green) line corresponds to (i_ and the lower one to a}+ within scheme 3. 
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4 Effective Dimer-Dimer Interaction 

Wc now turn our attention to the evaluation of an effective dimer-dimer interaction that takes both 
intra- and inter-layer interactions into account and uses the full dimer wave function as well. This 
procedure was discussed in |32j within the approximation scheme 2 outlined above. Here we elaborate 
and present a full comparison of various approximation and the influence on the effective dimer-dimer 
interaction. 



4.1 Derivation 

Consider two dimers in a bilayer. Let the coordinates of the polar molecules in each layer be ri, r2 in 
the first dimer and ra, r4 in the second dimer. Define the coordinates relative coordinates r — ri ~ r2 
and r' = — r4, and the center-of-mass (CM) for each dimer, R = (ri + r2)/2 and R' = {r^ + r4)/2. 
The distance between the CM of the two dimers is p = R — R! . The effective potential is obtained by 
integrating over the wave function of the dimers and over the CM coordinates with the condition that 
p = R R! (we specify the meaning of this below) . We have 



eff 



(p) = J drdr'dp'\^{ri,r2,r3,r4) 



[Vdip(ri - r2) + Vdip(r3 - r4) + V;iip(ri - r^) 
Vdip(r"i - Ti) + ydip(r2 - rs) + Vdip(r2 - r4)] , 



(18) 



where Vdip is the dipole potential in coordinate-space and W is the total wave function of the system. 
Here we have defined the total center-of-mass coordinates of all four molecules, p' = {R + R')/'2. For 
the latter we assume 



^{ri,r2,r3,ri) = ^{r)(b{r')^j{RWR'), 



(19) 



where (f) is the relative wave function and ip is the center-of-mass wave function of the dimers. We 
expect this to be a good approximation for a system of many particles when the dimer is strongly 
bound for large U and can be considered the effective constituent. We can write 



K//(p)= J drdr'dp'\<j>{r)\'\<p{r')\'^ 

mRrmR')\' 



V^{r) + V^ir') + V|| (p + '^^) 



V — V V + V T + T 

V\\{p-^^) + vAp + ^) + vAp-^ 



(20) 



where Vj_ is the inter-layer dipole potential and V\\ is the intra-layer dipole potential. Assume now 
that the center-of-mass part is proportional to 6{R R! p)6{{R + R')/'2 — Rcm), where Rqm is 
the center-of-mass position of the total four-body system which is unimportant here (we can fix the 
coordinate system so that Rcm = in any case). One can then drop the center-of-mass part of the 
integral. Likewise, the integral over the first two potential terms gives the potential energy of the dimer 
which is a constant that can be discarded. We are left with 



VM = / drdr'|0(r)|2|0(r')|' 



r r 

Vii{p+^—) 



r — r V + V 

+VAp-^—) + V^{p+^ 



V±{p- 



(21) 
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The Fourier transform is given by 

Kfr(fc) = JvMe^p-^'^Pdp 

drdr'\(l){r)\^\(j){r')\^ 21^1 (fe) cos ( fc 
-2yL(fc)cos (k ■ ^-—^ 



(22) 



where V||(/c) is the intra-layer and V±{k) is the inter-layer momentum-space potential respectively. 
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Fig. 6 (color online) Plot of the (normalized) effective dimer-dimer interaction, Kff(fc)/Veff(0), along with 
'F^{k) for [/ = 1, 2, 5, and 10, and with layer width w/d — 0.1. The lower of the two solid lines is the effective 

potential while the upper one is ^^{k) (and likewise for the dashed lines). The solid (blue) lines use the exact 
dimer wave function, whereas the dashed (black) lines are within the Gaussian approximation scheme 2. Note 
the differences in the horizontal scale of the panels. 

The explicit expression for the interaction potentials in momentum space are 

V\^iq)^^^{l^lF{\wq\)) (23) 
dV Zttw 

y^(q) = -27rZ?2|q|e-l9l^ (24) 

where F{x) = ^/^x[l - Erf(a;/\/2)]e^'/2 ^^^^ Erf(a;) the error function. We have assumed that only 
the ground-state of the optical lattice potential that creates the 2D geometry is occupied and we denote 
the transverse width by w. For the interlayer interaction we use the strict 2D hmit result which is a 
very good approximation for the values w < 0.3d to be used below. 
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4.2 Results 



From the expression in Eq. ( 22 ) , it is clear that to obtain the effective interaction of two dimers, 
one need the Fourier transform of the coordinate-space wave function squared. In Fig. [6j the effective 
interaction based on the exact wave function is compared to approximation scheme 2 for a layer width 
oi w = O.ld. As expected, the solid and dashed curves are clearly seen to approach each other as U 
increases. To explore the effect of the layer width, we plot Vcf{{k) in Fig. [T] for U = 5 and 10 for widths 
w — 0.2d and w ~ 0.3d. The agreement of the exact and approximate solution is striking except 
for U = 5 and w = 0.3d where the Gaussian wave function seems to somewhat underestimate the 
attraction at large kd. Note that for the physics to remain effectively 2D, we must have kw ^ 1 or 
kd ^ d/w. The range becomes smaller at larger w and at large kd the results are strongly modified 
by the 3D dipolar physics. 




We have demonstrated that the Gaussian approximation in scheme 2 is accurate at large U. In 
particular, this was the limit in which this approximation was used in |32j to calculate properties of 
the dimerized system at both zero and finite temperature. Furthermore, it is clear that this fact holds 
for a range of layer sizes, and the Gaussian approximation can be used to study the behavior of the 
roton instability of a bosonic dimer system that should emerge in the strong-coupling large U limit. 
The use of softer lattice that have larger w implies more attraction at large kd and the roton instability 
could become accessible at lower U values. 
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Fig. 8 (color online) Zero-temperature phase diagram with lines of vanishing chemical potential both with 
and without including the intralayer interaction as discussed in the text. Results based on the exact energy 
are shown as the thick soHd (blue) line. Also shown energies based on the the variationally optimized Gaussian 
(dashed (black) line), scheme 2 (solid (ref) line), and scheme 3 (dotted (green) lines). For comparison, the 
dash-dotted (purple) line to the right of the lines marked 'Including intralayer' shows the result of a cal cula tion 
of the crossover using the BCS gap equation including the self-energy term self-consistently (taken from [32]). It 
demonstrates the good agreement between BCS (traditionally weak-coupling theory) and the strong-coupling 
approach discussed in the current presentation. All calculations use w/d = 0.2. 



5 Many-body Physics 

We now consider the many-body physics of the bilayer at intermediate to strong coupUng where the 
BCS-BEC crossover is expected to happen. As discussed in the introduction, the two-component Fermi 
gas in 2D with zero-range interactions in the BCS approximation has a closed analytical solution |431 

mm 

A = ^/2EfEb (25) 
fi^EF~]^EB, (26) 

where A is the constant energy gap, fi the chemical potential (recall that is reserved for the reduced 
mass), Ep = h'^kp/2M the Fermi energy, and Eb the binding energy of the two-body bound state 
that the attractive zero-range potential supports. We define the crossover from weak-coupling BCS 
behavior to a two-body bound state BEC as the point where // = 0. If we write Eb = li^/MP, with 
I the size of the bound state, then /i = occurs when 2TrnP — 1, i.e. when the inter particle distance 
is roughly the size of the bound state. In our case with a bilayer of dipolar fermions, we expect these 
conclusions to be modified by the fact that we have long-range interactions that are both attractive 
and repulsive. 

From a Fermi liquid theory point of view, one could argue that since the long-range dipolar interac- 
tion is present on both side of the crossover, it provides a constant shift of the chemical potential and 
a modification of the effective mass entering the two-body bound that cancel each other and have no 
influence on the position where the crossover is expected. While this would be correct for small U, we 
are here considering the strong-coupling limit and this conclusion may no longer be true. In fact, recent 
studies have found that other many-body phases in the system such as the density-wave state are very 
strongly influence by the self-energy correction which drives the instability to larger values of U [301 
[35, 58, 60, 59, 61 > 62 . We expect qualitatively similar behaviour in relation to the pairing properties of 
the system. We therefore consider the problem from the strong-coupling point of view via an effective 
dimer approach. 

In a first approximation, we discard the intralayer interaction and consider the fl — Ep— ^|i?o| = 
line with Eq the binding energy of the bound state caused by the interlayer potential. Lines of vanishing 
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/2 can be seen to the right in Fig. |8] for the exact result and the different approximations to the binding 
energy. Except for scheme 3, the approximations are all very close the exact result, which is not 
surprising since the approximations are in one way or another derived from the demand that the size 
of the bound state is close to the exact result. 

The intralayer interaction must of course be taken into account, in particular in the strong-coupling 
regime. One way to proceed is to calculate the effective mass, m* and replace m — >■ m* in Ep as 
discussed above. This can be easily done in the weak-coupling regime where the result is analytic [34l 
|62,63'. However, we are interested in the strongly-coupled regime where we expect the dimer to be 
strongly bound and hence the relevant degree of freedom. We therefore consider instead the physics 
of an interacting Bose gas of two-body bound states for which the chemical potential to lowest order 
in many-body perturbation theory is jl — \nVca{^) |64j . The factor of one- half is needed since we are 
interested in the chemical potential per molecule and not per dimer. 

In order to include the intralayer interaction we consider now the chemical potential 

fL = EF + \7iV,n{Q)~]^\Eo\. (27) 

The lines where this expression vanishes are shown on the left side of Fig. |8] Clearly the extra repulsion 
at long distances (small kd) pushes the chemical potential up and for the binding energy to compensate 
this effect we need larger [/, i.e. the lines are pushed to the left to low densities. Again we see that 
the exact result and the approximations for the binding energy are very close. In Fig. [8] we have used 
w — 0.2d. Using a smaller w will move the lines even further left and vice versa for larger w. Tuning w 
could therefore be a necessary and convenient way to enlarge the BEC regime for experimental access. 

For comparison, we also show in Fig. [8] a calculation of the line of vanishing chemical potential 
using BCS theory, i.e. by solving the gap equations including both inter- and intralayer interactions 
as discussed in Ref. [32j. It is shown as a (purple) dash-dotted line on the right side of the strong- 
coupling approach including the intralayer term. The lines show that one finds agreement between the 
BCS approach and the strong-coupling approach discussed in the current paper. This suggests that 
the inclusion of the intralayer interaction through the ri,Vi,g(0) term in the strong-coupling expression 



of Eq. (27) is consistent with the BCS result. Note that the intralayer interaction is included self- 
consistently in the BCS equation with the full momentum-dependent potential [32) . The intralayer 
interaction indeed modifies the BCS phase and not merely by a shift of the chemical potential propor- 
tional to nV^ff(O). 

The use of lowest order perturbation theory to include the intralayer interaction (recall that the 
interlayer interaction is zero at long wavelength) corresponds to using the first Born approximation. 
This approximation traditionally fails in the strong-coupling regime. However, the crossover physics 
of interest happens at low density and we expect the first Born approximation to be reasonable. The 
perturbative parameter here is Ukpd which must be small. This is true in our phase diagram for 
?7 < 5. However, the interlayer interaction that is responsible for the bound-state is taken into account 



essentially exactly through the integration over the bound-state wave function in Eq. ( 22 ) . We therefore 
expect our results to still be at least qualitatively correct for larger U as well. 



5.1 Pauli Blocking Effect 

The above discussion was based on a binding energy for the dimer calculated from the pure two-body 
Schrodinger equation. The fact that the dimer must form in the proximity of other fermionic molecules 
was taken into account by adding the Fermi energy and the intralayer repulsion. We now proceed by 
a different route in order to the background Fermi sea into account when forming the dimer. This can 
be seen as a 'top-down' approach, where the Fermi sea is considered a background on top of which 
a dimer forms. In contrast, a 'bottom-up' approach would have to consider the full iV-body problem 
which is a much more difficult task. Our approach is the same as that used to demonstrate the binding 
of Cooper pairs in the case of short-range interactions for two particles in the presence of Fermi seas. 
Here it is known to be benificial to have the background, but this it is not a priory clear that the same 
holds with long-range dipolar interactions. Below, we calculate the energy in the presence of Fermi 
seas and show that the energy of the dimer goes up as the Fermi sea is increased. Thus taking the 
Pauli blocking effect into account to lowest order implies that the crossover could move in the phase 
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diagram. Notice that we have a Fermi sea background for both molecules that form the dimer. This is 
different from the polaron problem where only one of the particles has a background sea. Interetingly, 
the polaron problem in the bilayer setup with polar molecules has recently been studied by Klauwunn 
and Recati [55] . 

Consider the momentum-space two-body Schrodinger equation 

^'^^^^ + / J^^^^^ ' '^^'^(^^ ^ (^^^ 

where k is the relative momentum, which can be solved by discretization. We now include the effect 
of the background Fermi sea on the dimer by blocking all states with momenta less than the Fermi 
wave vector kp, i-e. we assume that 'F{k) = for |fc| < kp- We also assume that the dimer has zero 
center-of-mass momentum with respect to the Fermi sea, i.e. ki — — ^2 so that the relative momentum 
becomes the lab momentum. Here we are effectively assuming that the Fermi sea is inert in the sense 
that we neglect particle-hole pairs induced by the interactions. These correlations can be taken into 
account for instance through the variational ansatz employed for highly polarized Fermi gases with 
short-range interactions [66;, suitably modified to the case of a balanced system which is what we are 
concerned with here. This will be explored in future work. 
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Fig. 9 (color online) Dimer binding energy in units of /McF for [7 = 5, 10, and 15 as function of kpd. The 
solid (black) curve is the result of solving the momentum-space Schrodinger equation with all states below kp 
blocked. The dashed (blue) curve is the binding energy obtained by u sing a variationally optimized Gaussian 
without blocking plus twice the Fermi energy (according to Eq. (29l), while the dotted (red) curve includes 



twice the Fermi energy and also the lowest-order correction to the potential energy at finite kp, Eq. (32 1 



Before we discuss the numerical results, some limiting analytical expression can be obtained. Let 
us assume a Gaussian wave function which is a very good approximation for large U . Then the kinetic 
energy, T, in the Pauli blocked case is 

^^^[-<^'')^^ 

where h is the length scale of the Gaussian which can be obtained either variationally or within one 
of the approximation schemes above. Here we have taken care to re-normalize the wave function when 
the \k\ < kp is cut. This is an intuitively pleasing result, and we note that a similar expression can be 
obtained for the exponential ansatz. The kinetic energy per fermionic molecule can then be written 

2 7tf^-4M6^+^^- ^^^> 
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This should be related to the estimation of the chemical potential through /i = Ep — ^\Eq\ as shown on 
the right of Fig. |8] This approximation corresponds to using Pauli blocking in the kinetic energy while 
leaving the interaction energy term untouched. The matrix element of the potential in momentum 
space for the Gaussian wave function is 



(27r) 



-qd-b^q^ cos <, 



(31) 



where the factor e'-'^'^''-' comes from the re-normalization of the wave function, {^\^) = e '^''^'''> . In the 
kp ^ limit, the integral can be performed to yield the total energy given in Eq. ( 12 ). The corrections 



to the potential energy at finite kp can be obtained by expansion, which to lowest non-trivial order 
yields 



4 d3 



ikpb)' 



= - 


-2y 4 




V2n 


(^Cosh 


1 


+ Sinh 


1 











[V2y 



(32) 



(33) 



The function g{y) is positive so that the energy increases with increasing kp. 

In Fig. [9] we plot the exact solution including Pauli blocking (solid curve) for U = 5, 10, and 15 as 
a function of kpd. As one would expect from the arguments above, the energy increases with kpd in a 
quadratic fashion with a larger slope for larger U. For comparison, we have also plotted the results of 
a Gaussian variational calculation of the energy including the effect of blocking on the kinetic energy 
(dashed (blue) curve) and on both the kinetic and potential energy terms (dotted (red) curve). Recall 
again that the effect of blocking on the kinetic energy is exact, whereas that on the potential part 
is only to lowest non-trivial order in kpd. The dashed curve in Fig. [9j which includes the blocking 
in the kinetic energy, overestimates the binding. The inclusion of blocking in the potential (dotted 
curve in Fig|9]) is seen to provide an improvement at small kpd (particularly for larger U), but still 
underestimates the binding at large kpd. An expansion in large kpd can also be done but this is not 
analytical and provides only limited additional information. We will not pursue this further here. 

The fact that decrease of the binding energy with increasing size of the background Fermi sea is 
faster than predicted by simple addition of the Fermi energy implies that the phase diagram in Fig. [9] 
must be suitably modified to account for the background effects. Fig. [TO] presents a comparison of 
the lines of vanishing chemical potential when including the Pauli blocking through the exact solution 
of the momentum-space Schrodinger equation, both with and without including the effects of the 
intralayer interaction in the long-wavelength limit through addition of ^nVcs{0). This implies that the 
background environment reduces the importance of the dimer for the many-body physics. 



6 Conclusion and Outlook 

We have considered a bilayer system with fermionic polar molecules and in particular the two-body 
bound dimer state that the system supports at any dipole strength. The energy and wave function was 
calculated exactly and compared to various approximation schemes that provide convenient analytical 
expressions. We demonstrated that for dipole strength U > 1, a suitably chosen Gaussian provide a 
very good analytical approximation in both coordinate- and momentum-space. The effective dimer- 
dimer interaction was calculated within the different schemes and also found to be well approximated 
by using Gaussian two-body wave functions, which can be used, for instance, to study roton instabilities 
[32j . In conclusion, we find that the Gaussian choice is both an accurate and a convenient one due to 
its simple analytical properties. 

The corresponding many-body physics of the bilayer was studied in the context of the 2D Fermi 
gas where BCS-BEC crossover has been predicted for short-range interactions in the two-component 
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Fig. 10 (color online) Pliase diagram as in Fig.lSl The solid (black) lines obtained from the momentum-space 
Schrodinger equation including Pauli blocking of low momentum states with (right) and without (left) intralayer 
interaction, whereas the dashed (blue) lines are obtained by adding the Fermi energy to the unblocked binding 
energy. The latter corresponds to the solid (blue) lines in Fig. [s] Again we use w/d = 0.2. 



system. This crossover depends sensitively on the two-body bound state. The bilayer system considered 
in the present paper is very similar since the layer index can be mapped to a spin quantum number. 
By determining where the chemical potential of the dipolar system vanishes, we could map out a 
mean-field phase diagram of the BCS and (quasi)-BEC regions. We find this diagram to be largely 
independent of our approximation schemes for the binding energy. Here we took the effects of the 
background Fermi sea on the dimer into account both through a simple addition of the Fermi energy, 
but also by using Pauli blocking of all momenta below the Fermi momentum in the momentum-space 
Schrodinger equation for the dimer binding energy. This showed that blocking causes a decrease in the 
binding with the size of the Fermi sea that is faster than naively expected by addition of the Fermi 
energy to the non-blocked two-body binding energy. 

In the bilayer with fermionic polar molecules there are both repulsive and attractive parts. If we 
neglect the presence of the repulsive intralayer term in the bilayer, we recover the result from the short- 
range attractive 2D Fermi gas case. However, the long-range intralayer repulsive drastically changes 
the picture and pushes the BCS-BEC crossover to much smaller densities. This has clear analogs 
to electron- hole bilayers where the Coulomb interaction provides the long-range part. Here we have 
included the effect of Pauli blocking on the two-body bound-state energy and demonstrated that this 
pushes the crossover to even smaller densities, implicating that the BCS paired state occupies the 
majority of the mean-field phase diagram. 

To access the crossover region would require degenerate bilayer systems at very low densities. This 
is experimentally challenging since the critical temperature of this 2D setup is governed by the BKT 
transition [571l68j . The transition temperature is proportional to the superfluid density which can be 
depleted by both interaction and thermal effects. It is therefore driven to very small values at low 
density. In addition, we need dipole strengths in the range C/ ~ 1 — 10, which is higher than current 
experiments [H] . Using molecules with a larger moment like LiCs or a different optical lattice setup is 
therefore necessary to probe the crossover. Another possibility is for the system to develop a density 
wave 29,30,35 . However, strongly bound dimers should change this picture and it is not clear exactly 
where the density wave transition sits. In any case, we expect it to occupy the large-density and large- 
strength part of the phase diagram so the crossover physics should be accessible in the low-density and 
intermediate-strength region. 

Two-body states are one important constituent in bilayer system. However, 

due to the long-range character of the dipolar interaction, we can have three- or more-body states 
that are also bound, even though the particles are spatially separated and would not exhibit bound 
states with zero- or short-range interaction. This is a very interesting feature as it is very different from 
the paradigmatic crossover BCS-BEC from weakly bound Cooper pairs to strongly bound bosonic two- 
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body states which is driven by short-range two-body interactions. These questions wiU be pursued in 
future investigations. 
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